## Replication file 
library(foreign)
library(gtools)

setwd("/Replication files") #Set working directory where data is held

dfF <- read.dta("df_flanders_a.dta") # upload data for Flanders 


dfW <- read.dta("df_wallonia_a.dta") # upload data for Wallonia 

## Flanders
## change character vector to numeric values 
dfF$Evolution_economy <- as.numeric(dfF$Evolution_economy)
dfF$Patrimony_business <- as.numeric(dfF$Patrimony_business)
dfF$Patrimony_stocks <- as.numeric(dfF$Patrimony_stocks)
dfF$Patrimony_house <- as.numeric(dfF$Patrimony_house)
dfF$Patrimony_savings <- as.numeric(dfF$Patrimony_savings)
dfF$Immigrants_economy <- as.numeric(dfF$Immigrants_economy)
dfF$Immigrants_culture <- as.numeric(dfF$Immigrants_culture)
dfF$Immigrants_crime <- as.numeric(dfF$Immigrants_crime)
dfF$Marital <- as.numeric(dfF$Marital)
dfF$Education <- as.numeric(dfF$Education)
dfF$Employment <- as.numeric(dfF$Employment)
dfF$Party_close <- as.numeric(dfF$Party_close)
dfF$Interest_politics_general <- as.numeric(dfF$Interest_politics_general)

## Create variable for patrimony (high) and patrimony (low)
dfF$Patrimony_high <- dfF$Patrimony_business + dfF$Patrimony_stocks
dfF$Patrimony_low <- dfF$Patrimony_house + dfF$Patrimony_savings

## Create Age varibale 
dfF$Age <- 2019 - dfF$Birthyear

## Create index of trust variables
dfF$trust_index <- (dfF$Trust_courts + 
        dfF$Trust_police + 
        dfF$Trust_politicians + 
        dfF$Trust_polparties)/4

##Create index of immigraiton variables 
dfF$immigration_index <- (dfF$Immigrants_economy + dfF$Immigrants_culture + dfF$Immigrants_crime)/3


## Change the values of Government spendings for industry so higher values have positive meaning
dfF$Govspending_industry <- as.numeric(as.factor(dfF$Govspending_industry))
dfF$Govspending_industry1 <- ifelse(dfF$Govspending_industry > 5, NA, dfF$Govspending_industry)
dfF$Govspending_industry <- 6 - dfF$Govspending_industry1

## Wallonia
## change to numeric values 
dfW$Evolution_economy <- as.numeric(dfW$Evolution_economy)
dfW$Patrimony_business <- as.numeric(dfW$Patrimony_business)
dfW$Patrimony_stocks <- as.numeric(dfW$Patrimony_stocks)
dfW$Patrimony_house <- as.numeric(dfW$Patrimony_house)
dfW$Patrimony_savings <- as.numeric(dfW$Patrimony_savings)
dfW$Immigrants_economy <- as.numeric(dfW$Immigrants_economy)
dfW$Immigrants_culture <- as.numeric(dfW$Immigrants_culture)
dfW$Immigrants_crime <- as.numeric(dfW$Immigrants_crime)
dfW$Marital <- as.numeric(dfW$Marital)
dfW$Education <- as.numeric(dfW$Education)
dfW$Employment <- as.numeric(dfW$Employment)
dfW$Party_close <- as.numeric(dfW$Party_close)
dfW$Interest_politics_general <- as.numeric(dfW$Interest_politics_general)

## Create variable for patrimony (high) and patrimony (low)
dfW$Patrimony_high <- dfW$Patrimony_business + dfW$Patrimony_stocks
dfW$Patrimony_low <- dfW$Patrimony_house + dfW$Patrimony_savings

## Create Age varibale 
dfW$Age <- 2019 - dfW$Birthyear

## Create index of trust variables
dfW$trust_index <- (dfW$Trust_courts + 
        dfW$Trust_police + 
        dfW$Trust_politicians + 
        dfW$Trust_polparties)/4

##Create index of immigraiton variables 
dfW$immigration_index <- (dfW$Immigrants_economy + dfW$Immigrants_culture + dfW$Immigrants_crime)/3

## Change the values of Government spendings for industry so higher values have positive meaning
dfW$Govspending_industry <- as.numeric(as.factor(dfW$Govspending_industry))
dfW$Govspending_industry1 <- ifelse(dfW$Govspending_industry > 5, NA, dfW$Govspending_industry)
dfW$Govspending_industry <- 6 - dfW$Govspending_industry1

## Merge Flanders in Wallonia in one sample - Belgium
df <- rbind(dfF, dfW)

#####################
## Results manuscript
##################### 

#####################################################################################################
## Table 1. Ideological position of parties in Flanders and Wallonia based on respondents’ evaluation
#####################################################################################################

# Parties in Flanders (means and medians reported)

# Open VLD 
 summary(dfF$Leftright_1, na.rm=T) 
# N-VA
 summary(dfF$Leftright_2, na.rm=T) 
# Vlaams Belang
 summary(dfF$Leftright_3, na.rm=T) 
# CD&V
 summary(dfF$Leftright_4, na.rm=T) 
# PVDA
 summary(dfF$Leftright_5, na.rm=T) 
# Groen
 summary(dfF$Leftright_6, na.rm=T) 
# sp.a
 summary(dfF$Leftright_7, na.rm=T) 

# Parties in Wallonia (means and medians reported)

# Ecolo
 summary(dfW$Leftright_1, na.rm=T) 
# cdH
 summary(dfW$Leftright_2, na.rm=T) 
# MR
 summary(dfW$Leftright_3, na.rm=T) 
# PP
 summary(dfW$Leftright_4, na.rm=T) 
# DéFi
 summary(dfW$Leftright_5, na.rm=T) 
# PTB
 summary(dfW$Leftright_6, na.rm=T) 
# PS
 summary(dfW$Leftright_7, na.rm=T) 

#################
### Main Analysis 
#################

##############################################
## Multidimensional Economic Voting in Belgium
##############################################

## Regressions logit

model1 <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  									   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  									   Age + Sex + Marital + Education + Employment + Religion_frequency, data = df, family = binomial())

model2 <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = df, family = binomial())


#######################################################################
#### Predicted probabilities and marginal effect sizes based on model 1 
#######################################################################

library(margins)
marginal_effects <- margins(model1)
marginal_effects <- margins(model1, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

########################################################
## Figure 1. Multidimensional Economic Voting in Belgium
########################################################
library(ggplot2)

# Create a data frame with summary data based on summary(marginal_effects)
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0771, 0.0800, 0.0157, 0.0661),
  SE = c(0.0154, 0.0203, 0.0221, 0.0123),
  p_value = c(0.0000, 0.0001, 0.4784, 0.0000)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", "")))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
Belgium_AME <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.02, 0.105))

## Save plot in dedicated folder
ggsave("Belgium_AME.png", Belgium_AME, width = 13, height = 8, units = "in", dpi = 300)


############################
## Flemish elections of 2019
############################

###############################################
## Multidimensional Economic Voting in Flanders
###############################################

## Regressions logit

model1F <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())

model2F <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())

################################################################################
#### Predicted probabilities and marginal effect sizes based on model 1 Flanders
################################################################################

marginal_effects <- margins(model1F)
marginal_effects <- margins(model1F, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

##################################################################################################################
## Figure 2. Multidimensional Economic Voting across type of elections in Flanders - Federal Elections - left plot
##################################################################################################################

# Create a data frame with summary data based on summary(marginal_effects)
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0492, 0.0658, 0.0147, 0.0774),
  SE = c(0.0209, 0.0262, 0.0298, 0.0163),
  p_value = c(0.0182, 0.0119, 0.6233, 0.0000)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", "")))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
Belgium_AME_flanders <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.02, 0.105))

ggsave("Belgium_AME_flanders.png", Belgium_AME_flanders, width = 13, height = 8, units = "in", dpi = 300)


################################################################################
#### Predicted probabilities and marginal effect sizes based on model 2 Flanders
################################################################################

marginal_effects <- margins(model2F, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

####################################################################################################################
## Figure 2. Multidimensional Economic Voting across type of elections in Flanders - Regional Elections - right plot
####################################################################################################################

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0555, 0.0530, 0.0229, 0.0854),
  SE = c(0.0208, 0.0265, 0.0301, 0.0162),
  p_value = c(0.0076, 0.0453, 0.4462, 0.0000)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", "")))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
Belgium_AME_region_flanders <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) +
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.02, 0.105))

ggsave("Belgium_AME_region_flanders.png", Belgium_AME_region_flanders, width = 13, height = 8, units = "in", dpi = 300)


#############################
## 2019 Elections in Wallonia
#############################

###############################################
## Multidimensional Economic Voting in Wallonia
###############################################

## Regressions logit

model1W <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

model2W <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

################################################################################
#### Predicted probabilities and marginal effect sizes based on model 1 Wallonia
################################################################################

marginal_effects <- margins(model1W)
marginal_effects <- margins(model1W, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

##################################################################################################################
## Figure 3. Multidimensional Economic Voting across type of elections in Wallonia - Federal elections - left plot
##################################################################################################################

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0710, 0.0397, -0.0274, 0.0277),
  SE = c(0.0210, 0.0284, 0.0265, 0.0160),
  p_value = c(0.0007, 0.1627, 0.3002, 0.0831)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", "")))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
Belgium_AME_wallonia <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.06, 0.105))

ggsave("Belgium_AME_wallonia.png", Belgium_AME_wallonia, width = 13, height = 8, units = "in", dpi = 300)

##############################################################################
## Predicted probabilities and marginal effect sizes based on model 1 Wallonia
##############################################################################

marginal_effects <- margins(model2W, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

####################################################################################################################
## Figure 3. Multidimensional Economic Voting across type of elections in Wallonia - Regional elections - right plot
####################################################################################################################

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0627, 0.0129, 0.0153, 0.0425),
  SE = c(0.0225, 0.0345, 0.0297, 0.0188),
  p_value = c(0.0053, 0.7083, 0.6073, 0.0237)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", "")))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
Belgium_AME_region_wallonia <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) +
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.03, 0.105))

ggsave("Belgium_AME_region_wallonia.png", Belgium_AME_region_wallonia, width = 13, height = 8, units = "in", dpi = 300)

################################
### END of MAIN TEXT REPLICATION 
################################


##################
## Online Appendix
##################

#######################################
### Appendix 1. Full results Figure 1-3.
#######################################

#########################################################
##  eTable 1. Multidimensional Economic Voting in Belgium
#########################################################

model1 <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  									   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  									   Age + Sex + Marital + Education + Employment + Religion_frequency, data = df, family = binomial())

model2 <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = df, family = binomial())

#########################################################
## eTable 2. Multidimensional Economic Voting in Flanders
#########################################################

model1F <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())

model2F <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())


#########################################################
## eTable 3. Multidimensional Economic Voting in Wallonia
#########################################################

model1W <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

model2W <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

## End of Appendix 1.

#########################################
### Appendix 2. Linear Probability Models
#########################################

################################################################
##  eTable 2a. Multidimensional Economic Voting in Belgium (LPM)
################################################################

library(sandwich)
library(lmtest)

model1_lpm <- lm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										  Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										  Age + Sex + Marital + Education + Employment + Religion_frequency, data = df)

model2_lpm <- lm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										   Age + Sex + Marital + Education + Employment + Religion_frequency, data = df)

## Corrected standard Errors for the LPM models 
vv1 <- vcovHC(model1_lpm, type="HC1")
coeftest(model1_lpm, vcov = vv1)

vv2 <- vcovHC(model2_lpm, type="HC1")
coeftest(model2_lpm, vcov = vv2)

################################################################
## eTable 2b. Multidimensional Economic Voting in Flanders (LPM)
################################################################

model1F_lpm <- lm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										   Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF)

model2F_lpm <- lm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  											Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  											Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF)

## Corrected standard Errors for the LPM models 
vv1 <- vcovHC(model1F_lpm, type="HC1")
coeftest(model1F_lpm, vcov = vv1)

vv2 <- vcovHC(model2F_lpm, type="HC1")
coeftest(model2F_lpm, vcov = vv2)

################################################################
## eTable 2c. Multidimensional Economic Voting in Wallonia (LPM)
################################################################

model1W_lpm <- lm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										   Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW)

model2W_lpm <- lm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  											Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  											Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW)

vv1 <- vcovHC(model1W_lpm, type="HC1")
coeftest(model1W_lpm, vcov = vv1)

vv2 <- vcovHC(model2W_lpm, type="HC1")
coeftest(model2W_lpm, vcov = vv2)

## End of Appendix 2.

###################################################
## Appendix 3. Restricted Models on the Sample Size
###################################################

#############################################################################################
## eTable 3a. Multidimensional Economic Voting in Flanders (Restricted models on sample size)
#############################################################################################

### Resttricted Flanders sample  
library(dplyr) 

# select variables from regression
dfFR <- dfF %>% select(Vote_federal_incumbent,
                     Vote_regional_incumbent,
                     Evolution_economy,
                     Patrimony_high,
                     Patrimony_low,
                     Income_government_action,
                     trust_index,
                     Europeanisation,
                     immigration_index,
                     Party_close,
                     Leftright_self,
                     Interest_politics_general,
                     Satisfaction_fedgov,
                     Age,
                     Sex,
                     Marital,
                     Education,
                     Employment,
                     Religion_frequency)
dfFR_NA <- na.omit(dfFR) # omit all na's 

# regression models on restricted sample in Flanders
regFnR <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + trust_index + Europeanisation + immigration_index + 
  									   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  									   Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfFR_NA, family = binomial())

regFrR <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfFR_NA, family = binomial())




#############################################################################################
## eTable 3b. Multidimensional Economic Voting in Wallonia (Restricted models on sample size)
#############################################################################################

### Resttricted Wallonia sample  

# select variables from regression
dfWR <- dfW %>% select(Vote_federal_incumbent,
                     Vote_regional_incumbent,
                     Evolution_economy,
                     Patrimony_high,
                     Patrimony_low,
                     Income_government_action,
                     trust_index,
                     Europeanisation,
                     immigration_index,
                     Party_close,
                     Leftright_self,
                     Interest_politics_general,
                     Satisfaction_fedgov,
                     Age,
                     Sex,
                     Marital,
                     Education,
                     Employment,
                     Religion_frequency)
dfWR_NA <- na.omit(dfWR) # omit all na's 

# regression models on restricted sample in Wallonia
regWnR <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + trust_index + Europeanisation + immigration_index + 
  									   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  									   Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfWR_NA, family = binomial())

regWrR <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfWR_NA, family = binomial())


#######################################################################################
## eFigure 2. Predicted Probability of Multidimensional Economics by item and Vote for
## incumbent government in 2019 elections in Belgium by federal and regional elections.
#######################################################################################

################################################################################################################################################
## Average Marginal Effedcts by value of  Economic evaluation (valence), patrimony and Income government action (positional) - Federal Elections
################################################################################################################################################

library(effects)

## Extract effects in object for ploting based on model1 
marginal_effects_Evolution_economy <- effect("Evolution_economy", model1)
marginal_effects_Patrimony_high <- effect("Patrimony_high", model1)
marginal_effects_Patrimony_low <- effect("Patrimony_low", model1)
marginal_effects_Income_government_action <- effect("Income_government_action", model1)

#################################################################
## Economic evaluation - Federal elections - top plot on the left 
#################################################################
library(ggplot2)

plot_evo <- ggplot(data = as.data.frame(marginal_effects_Evolution_economy), aes(Evolution_economy, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Economic evaluation", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Evolution_economy)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_valence <- plot_evo + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_valence.png", Belgium_valence, width = 8, height = 6, units = "in", dpi = 300)

#####################################################
## Patrimony (high) - Federal elections - left middle
#####################################################

plot_pat_high <- ggplot(data = as.data.frame(marginal_effects_Patrimony_high), aes(Patrimony_high, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Patrimony (high)", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Patrimony_high)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_patrimony_high <- plot_pat_high + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_patrimony_high.png", Belgium_patrimony_high, width = 8, height = 6, units = "in", dpi = 300)

####################################################
## Patrimony (low) - Federal elections - left middle
####################################################

plot_pat_low <- ggplot(data = as.data.frame(marginal_effects_Patrimony_low), aes(Patrimony_low, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Patrimony (low)", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Patrimony_low)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_patrimony_low <- plot_pat_low + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_patrimony_low.png", Belgium_patrimony_low, width = 8, height = 6, units = "in", dpi = 300)

#############################################################
## Income government action - Federal elections - left bottom   
#############################################################

plot_position <- ggplot(data = as.data.frame(marginal_effects_Income_government_action), aes(Income_government_action, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Income government action", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Income_government_action)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_positional <- plot_position + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_positional.png", Belgium_positional, width = 8, height = 6, units = "in", dpi = 300)

################################################################################################################################################
## Average Marginal Effedcts by value of  Economic evaluation (valence), patrimony and Income government action (positional) - Regional Elections
################################################################################################################################################

## Extract effects in object for ploting based on model2 
marginal_effects_Evolution_economy <- effect("Evolution_economy", model2)
marginal_effects_Patrimony_high <- effect("Patrimony_high", model2)
marginal_effects_Patrimony_low <- effect("Patrimony_low", model2)
marginal_effects_Income_government_action <- effect("Income_government_action", model2)

###################################################################
## Economic evaluation - Regional elections - top plot on the right 
###################################################################

plot_evo <- ggplot(data = as.data.frame(marginal_effects_Evolution_economy), aes(Evolution_economy, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Economic evaluation", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Evolution_economy)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_valence_region <- plot_evo + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_valence_region.png", Belgium_valence_region, width = 8, height = 6, units = "in", dpi = 300)

#######################################################################################
## Patrimony (high) - Reginal elections - right middle
#######################################################################################

plot_pat_high <- ggplot(data = as.data.frame(marginal_effects_Patrimony_high), aes(Patrimony_high, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Patrimony (high)", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Patrimony_high)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_patrimony_region_high <- plot_pat_high + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_patrimony_region_high.png", Belgium_patrimony_region_high, width = 8, height = 6, units = "in", dpi = 300)

######################################################################################
## Patrimony (low) - Reginal elections - right middle
#######################################################################################

plot_pat_low <- ggplot(data = as.data.frame(marginal_effects_Patrimony_low), aes(Patrimony_low, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Patrimony (low)", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Patrimony_low)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_patrimony_region_low <- plot_pat_low + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_patrimony_region_low.png", Belgium_patrimony_region_low, width = 8, height = 6, units = "in", dpi = 300)


#######################################################################################
## Income government action - Regional elections - left bottom
#######################################################################################

plot_position <- ggplot(data = as.data.frame(marginal_effects_Income_government_action), aes(Income_government_action, fit)) +
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "gray") +
  labs(x = "Income government action", y = "", title = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 26, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 26, family = "Times New Roman", face = "bold"))

predicted_values <- as.data.frame(marginal_effects_Income_government_action)
predicted_values$label <- round(predicted_values$fit, 2)

Belgium_positional_region <- plot_position + geom_text(data = predicted_values, aes(label = label), size = 8, vjust = 2) +
  scale_y_continuous(limits = c(0, 1))

ggsave("Belgium_positional_region.png", Belgium_positional_region, width = 8, height = 6, units = "in", dpi = 300)

## End of Appendix 3

#########################################################
## Appendix 4. Alternative dependent variable in Wallonia
#########################################################

# Recode variable to center-right/right parties in Wallonia
dfW$Vote_federal_right <- ifelse(dfW$Vote_federal == 3, 1, ifelse(dfW$Vote_federal == 2, 1, ifelse(dfW$Vote_federal == 4, 1, ifelse(dfW$Vote_federal == 5, 1, ifelse(dfW$Vote_federal == 8, NA, 0)))))
dfW$Vote_regional_right <- ifelse(dfW$Vote_regional == 3, 1, ifelse(dfW$Vote_regional == 2, 1,  ifelse(dfW$Vote_regional == 4, 1, ifelse(dfW$Vote_regional == 5, 1, ifelse(dfW$Vote_regional == 8, NA, 0))))) 

#############################################################
## eTable 4a. Vote for center-right/right parties in Wallonia
#############################################################

model1W_right <- glm(Vote_federal_right ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										  Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										  Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

model2W_right <- glm(Vote_regional_right ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										   Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										   Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

###################################################################
## eTable 4b. Vote for center-right/right parties in Wallonia (LPM)
###################################################################

model1W_right_lpm <- lm(Vote_federal_right ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  											 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  											 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW)

model2W_right_lpm <- lm(Vote_regional_right ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  											  Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  											  Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW)

vv1 <- vcovHC(model1W_right_lpm, type="HC1")
coeftest(model1W_right_lpm, vcov = vv1)

vv2 <- vcovHC(model2W_right_lpm, type="HC1")
coeftest(model2W_right_lpm, vcov = vv2)

## End of Appendix 4

#####################################################################################
#####################################################################################
## Appendix 5. Out-of-sample test of hypotheses. -- Please see Replication_code_2.txt
#####################################################################################
#####################################################################################

##########################################################################################
## Appendix 6. Government spending for companies and industry across Flanders and Wallonia
##########################################################################################

###########################################################
## eTable 6. Government spending for companies and industry
###########################################################

# linear regression
reg_ind <- lm(Govspending_industry ~ as.factor(Patrimony_high) * as.factor(Region) + as.factor(Patrimony_low) * as.factor(Region) + Evolution_economy + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
									 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
									 Age + Sex + Marital + Education + Employment + Religion_frequency, data = df)

summary(reg_ind)

## End of Appendix 6

############################################################################################
## Appendix 7. Average Marginal Effects for Figure 2 and 3 based on full model specification
############################################################################################

#####################################################################################
## eTable 7a. Average Marginal Effects calculated based on Flanders Federal elections
#####################################################################################

model1F <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())

marginal_effects_F <- margins(model1F)
summary(marginal_effects_F)

#####################################################################################
## eTable 7b. Average Marginal Effects calculated based on Wallonia Federal elections
#####################################################################################

model1W <- glm(Vote_federal_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_fedparl + trust_index + Europeanisation + immigration_index + 
  										Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())

marginal_effects_W <- margins(model1W)
summary(marginal_effects_W)

######################################################################################
## eTable 7c. Average Marginal Effects calculated based on Flanders Regional elections
######################################################################################

model2F <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfF, family = binomial())

marginal_effects_F2 <- margins(model2F)
summary(marginal_effects_F2)

#####################################################################################
## eTable 7b. Average Marginal Effects calculated based on Wallonia Federal elections
#####################################################################################

model2W <- glm(Vote_regional_incumbent ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Trust_reggov + trust_index + Europeanisation + immigration_index + 
  										 Party_close + Leftright_self + Interest_politics_general + Satisfaction_fedgov + 
  										 Age + Sex + Marital + Education + Employment + Religion_frequency, data = dfW, family = binomial())


marginal_effects_W2 <- margins(model2W)
summary(marginal_effects_W2)

## End of Appendix 7.


### End of replication file.